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Abstract 

We investigate the effects of higher order curvature corrections to Einstein's Gravity on the 
critical phenomenon near the black hole threshold, namely the Choptuik phenomenon. We simulate 
numerically a five dimensional spherically symmetric gravitational collapse of massless scalar field 
in Einstein-Gauss-Bonnet (EGB) gravity towards a black hole formation threshold. When the 
curvature is sufficiently large the additional higher order terms, affect the evolution of the whole 
system. Since high curvature characterizes the region when the critical behavior takes place this 
critical behavior is destroyed. Both the self similarity and the mass scaling relation disappear. 
Instead we find a different behavior near the black hole threshold, which depends on the coupling 
constant of the higher order terms. The new features include a change of the sign of the Ricci 
scalar on the origin which indicates changes in the local geometry of space-time, and never occurs in 
classical general relativity collapse, and oscillations with a constant rather than with a diminishing 
length scale. 
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I. INTRODUCTION 



Gravity is described by Einstein's theory of general relativity (GR), which predicts the 
existence of Black holes - trapped regions from which nothing can escape. Black holes can 
form from regular initial data that do not contain a black hole already. Isolated system 
in GR can end up in two qualitatively different states. Data that forms a black hole in 
the evolution and data that doesn't and in which the mass-energy disperses to infinity. A 
simplest possible model for such a system is the collapse of a spherically symmetric minimally 
coupled massless scalar field. The final fate - whether it collapses or not - depends on the 
"strength" of the initial data. 

In a pioneering work Choptuik [1] explored the transition between the two regimes. He 
discovered that the black hole threshold shows both surprising structure and surprising 
simplicity. Universality, power-law scaling of the black hole mass, and scale echoing have 
given rise to the term "critical phenomena". [2, 3]. Choptuik [1] considered a one parameter 
families of initial data describing a collapsing scalar field (see figure ?? for the structure 
of space time)/ He have shown that for each family of initial data parametrized by p (for 
example amplitude of the initial pulse) there exist a critical value p*. For p > p* we have a 
supercritical collapse and a black hole forms. For p < p* the collapse is subcritical and the 
field disperses to infinity leaving a flat space. Choptuik gave a highly convincing numerical 
evidence that by fine-tuning the parameter p to the threshold value p* an infinitely small 
black hole can be created. 

The critical solution itself is universal. For a finite time in a finite region of space the 
spacetime converges to one and the same solution independent of the initial data. The 
critical solution is discretely self similar (DSS), namely, it is invariant under rescaling by a 
particular finite factor, or its integer powers. Let Z*(r, t) be the critical solution (collectively 
for all the parameters - the scalar field and the metric), the critical solution is the same when 
rescaling space and time by factor e A : 

Z*(r,t) = Z*{re A ,te A ). (1) 

The field and metric functions pulsate periodically with ever decreasing temporal and spatial 
scales, until a black hole forms in supercritical collapse, or the field disperses in subcritical 
collapse (see figure 2). This universal phase ends when the evolution diverges towards a black 
hole formation or towards dispersion, depending on whether p > p* or not. In supercritical 
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collapse above p* arbitrary small black holes are formed as p — > and the black hole mass 
scales as a power law: M oc (p — p*) 7 , where 7 is universal. It depends on the type of 
collapsing matter (and on the dimension) but it is independent of the initial data family. 
Similar critical phenomena were found in many other types of matter coupled to gravity, 
with spherical symmetry and beyond it (See e.g. [2] for a review). The echoing period A and 
critical exponent 7 depend on the type of matter and on the dimension, but the phenomena 
appears to be generic. 

Einstein's equations are derived from the Hilbert action, which is linear in the Ricci 
scalar - R. It is natural to expect that higher terms in R will appear in a more general 
theory and it is interesting to explore their possible role. To do so we have to explore a high 
curvature regions of space time where such terms are significant. Black hole formation is 
a natural place to do so as spacetime becomes highly curved as the matter fields collapse. 
This behavior usually takes place near the singularity, which is typically hidden inside the 
black hole. However they also appear near the threshold for black hole formation when the 
Choptuik phenomenon take place. Therefore we explore here the gravitational collapse of 
a spherically symmetric massless scalar field with higher order corrections to the Hilbert 
action. 

Addition of even the simplest R 2 term induces 4 th order derivatives of the metric in 
the resulting equations of motion. The original Einstein equations are second order in the 
metric and the additional generic 4 th order corrections will govern the equations and change 
completely their character. To overcome this we consider a special case. Lanczos [4, 5] 
found a generalization of Hilbert's Lagrangian which is quadratic in the Riemann tensor 
and its contractions, but its variation yields a system of equations that, like Einstein's, is 
second order in the metric derivatives. This correction to the Lagrangian, called the Lanczos 
Lagrangian or the Gauss-Bonnet term [6, 7], is given by: 

L GB = R 2 - AR^R"" + R^R^. (2) 

The combined theory that includes the Hilbert Lagrangian and the Gauss-Bonnet term is 
called Einstein-Gauss-Bonnet (EGB) gravity. In four dimensions the Gauss-Bonnet term is a 
pure divergence, just like Hilbert's Lagrangian in two dimensions, and it does not contribute 
to the field equations in four dimensions. To overcome this we consider here gravitational 
collapse in 5 dimensional space-time, which is the simplest system in which the Gauss-Bonnet 
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FIG. 1: The Penrose diagram in GR of the space time that is expected to form in a gravitational collapse 

of a shell of in-falling scalar field to a black hole. When the shell is far away from the origin the self 
gravitational effects are small. When it comes closer to the origin gravitational field becomes stronger. If 
the field doesn't collapses to a black hole the diagram remains Minkowsky- flat, and the event horizon or 
the singularity don't exist, of course. Light blue lines indicate the numerical domain of integration used in 
the current work. The initial hypersurface is a null ray. Since the field is massless it behaves light-like and 

it moves along null rays. 

term contributions can affect the evolution [8] . 

We present here the results of a numerical investigation of the influence of higher order 
curvature correction, namely the Gauss-Bonnet term, on the properties of spherically sym- 
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(a) Contours of the scalar field function svsw and v (b) The field function s at the origin (r = 0) vs u. 

FIG. 2: The field function s in classical GR in 5 dimensions, for a slightly subcritical collapse. The field 
oscillates with DSS pattern. The field pulsations decreasing in temporal and spatial scales. 

metric scalar field collapse in 5 dimensions. In particular we focus on the behavior of the 
Choptuik critical phenomenon. The structure of the paper is as follows: in §11 we describe 
the overall model and the basic equations. In §111 we discuss the numerical structure and the 
numerical difficulties that arise in the calculations. Simulation results and their discussion 
are presented in §IV. 

II. THE MODEL AND THE BASIC EQUATIONS 

We consider the collapse of a spherically symmetric massless scalar field, in 5 dimen- 
sional space time described by the metric g a p. We use units in which G — c — 1. 

A. The Metric 

The scalar field is massless and it propagates along the light cone. Hence we describe the 
5 dimensional asymptotically flat spacetime in double null coordinates: 

ds 2 = —a 2 (u, v)dudv + r 2 (w, v )dQl, (3) 
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where our coordinates are: (w, t>, y9, 6>, 7) and d£l\ = [sin 2 ^ (d6 2 + sin 2 9d(p 2 ) + d^y 2 ] is the 
metric on 3 dimensional unit sphere (see [9, 10] for an alternative Hamiltonain formulation 
of this problem). The coordinate u is the retarded time coordinate and a constant u describes 
an outgoing null trajectory. Similarly v is the advanced time coordinate and surfaces with 
a constant v are the ingoing null trajectories, r = r(u,v) is the area coordinate and r = 
is the origin of the spherical symmetry. This definition of the metric is unique only up to a 
change of variables of the form v — > v(v), u — » u(u). This gauge freedom will be fixed later 
by the choice of the initial conditions. 

B. The Scalar Field 

The Lagrangian density of the field is: 



L = 



(4) 



and the corresponding equation of motion is: 



; V = o. 



(5) 



The energy-momentum tensor is given by: 
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(6) 



C. The equations of Motion 

The overall action that includes the Hilbert action and the Gauss-Bonnet term is: 

1 



J d n x^~g 



16tt 



(R + aL GB ) 



+ S, 



matter i 



(7) 



where a is the coupling constant and it has dimensions of (length) 2 . The value of the 
coupling constant a is unknown, we assume that < a « 1. As we are interested in a 
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values which are significant and influence the solution. Therefore we will look for values 
that are large enough so that the correction term which is of order aR 2 is comparable or 
larger then the Hilbert term R. Namely we will be interested in cases where aR > 1. 
The corresponding gravitational field equations are: 

+ aH^ = n 2 n T^ u , (8) 



where G^ v is Einstein tensor: 



G>„ — Rfj, u — -Rg^u, (9) 



T M1 , is the energy-momentum tensor given in (6), and 

H^ u = 2 [RR^v — 1R m R a v — 2R a/3 R^ avl3 + R^ 1 R va ^\ — ^Q^Lgb- (10) 

We convert the equations to a set of first order differential equations. To do so we define: 

s = V4^G(f) (11) 
z = s v (12) 
w = s u (13) 

d EE ^ (14) 

(X 

f = r u (15) 

9 = r v , (16) 

where = J^. We obtain, using equation (8) four independent first order equations: 

f = 9 f 

a 3 [a 2 + 4ar7] 

2GraV 
3 [a 2 + 4ar?] 

a 2 rr] 
'2 [a 2 + 4m?] 

and 



a n 2Gra?w 2 

/„ = 2/ — f ; mm component; (17) 



2GraV 

^ = 2ga — - r 2 — - — y ; component; (18) 
a 2 rr] 

fv = 9u = - r g i : — t ! ^ component; (19) 



, -9a 2 r/[(4ar/) 2 -3a 4 ]+4a 2 G^[-9(a 2 + 4ar/) 2 + 32a 2 G^d a 2 

<*» = 36(a 2 + 4^)3 '^77 + ^ components. 

(20) 
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The scalar field equation derived from (5) is: 



W r . 



(21) 



We have defined here an auxiliary function 77: 

a 2 + 4fg 

V = — -0 — • 



(22) 



This is useful to stabilize the numerical solution, as explained later in Section III. The 
evolution equation for r\ is: 

2r]{g-rd) 2a 2 (3gr] + Afz 2 ) 



Vv = 



(23) 



r 3r(a 2 + Aarf) 

Equations - (17), (18), (19), (20), (21), (23) - together with (12), (13), (14), (15) and (16) 
are a set of eleven coupled first order differential equations that describe the system. Since 
there are only nine functions - r, a, s, f, g, z, w, d and r\ - two equations are redundant and 
are simply the usual constraint equations. We choose those to be Eqs. (15) and (17). They 
are not evolved in the integration, but they are monitored to verify that they are indeed 
satisfied during the evolution. 

The original GR equations are invariant under rescaling. However, the coupling constant 
a has units of length squared [L 2 ]. It's introduction to the system destroys this scale 
invariance. Under the rescaling: r — > pr the equations behave as: 
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(24) 

(25) 
(26) 



d u 



-9a 2 r][{4ar]) 2 - 3a 4 ] + 4a 2 Gzw[-9{a 2 + Aar]) 2 + 32a 2 Gzwa] 



36(a 2 + 4ar]) 3 

1 -9a 2 V [(^) 2 - 3a 4 ] + 4a 2 Om>[-9(a 2 + + 32a2 ^ a ] 



36(a 2 + 4 -^f 



(27) 



Obviously the scale invariance is broken. At smaller scales (p — > 0) the non-invariant 
elements, which rescale as p~ 2 , become dominant and govern the equations. This leads to 
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the deviations from the classical GRB behaviour that we demonstrate numerically later. In 
particular the deviations from the GR behavior happen near the critical point where the 
curvature is large and the non scale invariant Gauss-Bonnet terms are large. These terms 
eventually destroy the classical Choptuik phenomenon both when black holes form and when 
they don't. This happens independently of the value of a. 



D. The Boundary Conditions 

Regularity and differentiability at the origin r = require the following boundary condi- 
tions: 

g = -f = l d r s = 0, d r a = 0. (28) 

These conditions imply 

a v = a u ; w = z ; r] = (—3a + v^a 2 + A&az 2 ) (29) 

12a 

at the origin. 

A second set of boundary conditions is set implicitly at infinity. This is trivial in GR 
where a Schwartzchild has naturally an asymptotically flat space time. However, [11] have 
shown that a Schwartzchild black hole in EGB gravity can have either asymptotically flat 
or AdS structure. Since we consider only a finite region of space time (see figure 1) we don't 
examine here to which to the two branches the collapsing black hole will lead. 



E. The Ricci Scalar 



The Ricci scalar curvature, R, is given by: 

-9a 6 wz + 432?7 4 a 3 + 72a 2 r) 2 (r) - 2wz)a 2 + a 4 (32wV - 72wzrj - 27rf)a 
R ~ 8 902(02 + 4^)3 ' (30J 

The Ricci curvature describes the local geometry of the space-time. The value of the Ricci 
scalar at the origin is of special interest: 

R(r = 0) = -. 16Z " + - I - I I . (31 ) 

0^902 + 48^0 a y^ a 2 + iA z 2 a 



8z 2 

If a = the Ricci scalar at the center is always negative: For q^O the situation is 



a 2 



more complicated. The Ricci scalar is negative while the following condition is satisfied: 

45a 2 

« < iSp- (32) 

The simulations show that along the evolution of the collapsing field system, the metric 
function - a(u,v) approaches zero and z(u,v), the field derivative, is growing to a very big 
values, the closer we are to the critical amplitude, the larger is the value that z approaches. 
Therefore, the condition of eq.(32) is violated at some point and the Ricci scalar changes 
sign. This change of sign heralds the deviation from the classical behavior. 



F. The Black Hole Mass 

The analysis of mass scaling relation in the critical phenomenon requires a function for 
a black hole mass. The ADM mass of a black hole in higher dimensional GR [12] is: 

when r s is Schwarzschild radius, D - the dimension, Gr> is the D-dimensional Newton con- 

„ O-l 

stant and A D _ 2 is the area of a unit sphere: A D _ 2 = . D-1 . ■ However, in EGB there is an 

M 2 ) 

additional term and for D = 5 dimensions the ADM mass is [11, 13-15] : 



37r r,. 2a, 

M =8G r ' (1+ # (34) 



The second term in this equation implies that as the size of the black hole decreases (r s — > 0) 
its mass in EGB approaches a constant positive value, M — > M = 3ira/ 'AG [14, 15]. This 
implies that there is a mass gap and all black holes (for a > 0) must have an ADM mass 
larger than M . One can resort to a different definition of the black hole's mass and instead 
of using the ADM mass one can calculate the mass of the apparent (trapping) horizon using 
an EGB quasi-local mass [16, 17]. Avoiding this problem we will consider, for simplicity, in 
the following the scaling of the black hole's radius instead of the scaling of the black hole's 
mass. 

We define the critical exponent 7 such that \p — p*| 7 has a dimension of length. Instead 
of examining the dependence of the ADM mass on p we will examine the dependence of the 
black hole's radius, r s . We expect following the GR case to find r s oc (p — p*) 1 . 
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G. Initial conditions 



We turn to discuss the initial conditions. We consider here the gravitational collapse 
of a shell of in-falling scalar field. Figure 1 shows the domain of the current numerical 
work embedded inside an expected Penrose diagram. The u-v plane is covered by a two 
dimensional grid, as described in Figure 3. The origin r = is included in the domain, and 
it is chosen to be at u = v. Therefore the relevant part of u-v space is v > u. Our metric is 
defined up to a coordinate transformation, as it was mentioned in section II A. This gauge 
freedom is fixed by specifying the metric functions on the initial hypersurface, an initial ray 
with a constant retarded time u = Ui = 0. In flat - Minkowsky - space-time the conventional 
definition of the null coordinates is u = t — r and v = t + r. Since we are dealing with 
spherical shell, the space-time is flat in two regions - inside the shell and at asymptotically 
large radii. The integration starts far away from the event horizon and therefore the metric 

is nearly flat. Thus we define the area coordinate r along the initial null surface u — Ui — 

v 

as in a flat space time r = -. The fact that the spacetime is only approximately Miknowski 
is pronounced by the deviation of other metric function - a(u,v) - from its flat-space value 
(a = 1). In addition, we set a = 1 at the origin, r = 0, at one point u — v — 0. From 
here we can obtain all the other functions on the initial hypersurface by integration from 
the origin. 

The exact shape of the initial scalar field is unimportant as the Choptuik behavior is 
universal and independent of this shape. We choose the initial scalar field profile along an 
outgoing hypersurface u = Ui = to be a Gaussian: 



where the constants v c and a determine the initial position and width of the shell and the 
constant p is the amplitude of the pulse, p is the strength parameter of the initial data and 
it is the dynamical parameter that we vary to explore the Choptuik phenomenon. 

Having specified the functions r and s on the initial hypersurface, we can derive analyt- 
ically z and g, which are simply the derivatives z = s v and g — r v — 1/2. All the other 
functions on the initial hypersurface - /, a, w, d and r\ - are obtained by integrating the 
appropriate equations from the origin. 




(35) 
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III. NUMERICAL METHODS 



We follow here the methods developed by Sorkin and Oren [3]. However, as discussed 
later, further steps, including the addition of the variable i] are needed here to stabilize the 
code near the origin. The problem arises, because of the stronger non linear behavior of the 
Gauss-Bonnet terms. 

A. The Integration Sceme 

Our domain of integration is a equilateral right angle triangle in a u-v plane: < u < 
u m ax, < v < v max and v > u. v max = u max is chosen to cover the interesting relevant region. 
An illustration of the domain is sketched in figures 1 and 3. The simplest computational 
cell is square with grid spacing h u = h v = h. Triangular cells near the origin {u = v) are 
treated separately. 

The integration begins from the lower line of constant u = and propagates to the next 
line. Once the solution along an outgoing hypersurface with constant u value, u = U — h, 
is known, d and z are propagated to the next line, u = U, using equations (20) and (21) 
correspondingly. Then equations (21), (19), (18), (12), (16), (14) and (23) are integrated 
using a fourth order Runge-Kutta algorithm from the origin outward along v to obtain the 
functions w ,f, g, s, r, a and i], respectively. The remaining equations are not used directly, 
but they must be satisfied and are used to test the numerical solution. 

The first points near the origin are treated separately. The region near the origin is 
unstable. The instability arises near the origin where discretization errors are amplified, 
especially in source terms that involve a division by r. To resolve this problem we take 

several steps. First we introduce the new variable r\ = . This, unnecessary from 

the first glance, variable is an algebraic combination of others. It appears in every source 
function. It includes a division by r 2 which is very sensitive to errors near the origin. The 
independent evolution of r\ using equation (23) help stabilizing the source functions. In 
addition we have to take a few more steps: We use at points near the origin a more stable, 
second order Runge-Kutta algorithm. 

• Instead of integrating the function functions / and w along v we evaluate these func- 
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U tUjo,. ,V Vj- 



u=0,v=0 



u=0,v=v 



FIG. 3: The domain of integration. The calculations employs the two previous lines L\ and L2 in 
addition to the line, L 0l that is currently being solved. The boundary conditions on r = involving d r are 
implemented using the 3-point derivatives along the diagonal line. Smoothing of some functions near r = 
(at a point marked by cross) is done using past light cone points (marked by circles). Figure taken from [3]. 



tions using a Taylor expansion, e.g.: 



{dvf 



f(v) = f(v ) + dvf v (v ) + K -^fw(v ) + 0(dv) 6 , 



when Vq is the v value on the origin and dv = v — v$. 



(36) 



Additionally we smooth the functions z and d. First we evaluate the function, at some 
point P, then its value is smoothed with the values of the same function on points on 
the past light cone of P (see Figure 3). For example, for the function z at the point 
marked by cross in Figure 3 the new, smoothed, z is calculated according to: 



u- z, 



1 



3 + w' 



(37) 



where z e is the value obtained from the evolution equation, Z{ are the extrapolated 
values of z along the 3 directions of the past light cone: Zi(u) = 2z{u — h) — z(u — 2h), 
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and co is a weight parameter, which is varying for different functions and different code 
parameters, but typically u oc 0.1. 

Every one of those actions on its own stabilizes the integration but is insufficient to keep 
the code completely stable till the end. For the classical GR evolution, with a = 0, rj is not 
needed, however it is essential for the more general evolution. This combination prevents 
the code from crashing at least for low enough values of a. However, a strong penalty is 
paid as the combination and in particular the smoothing reduces the convergence rate of 
the code to a linear order. Using this algorithms we are able to get a stable and convergent 
evolution for small values of a. However, the code still becomes unstable for large values of 
a, usually when the field amplitudes that are close to the critical one. 



B. Numerical tests 

We performed a series of simulations with step sizes h, \ and \ in order to determine the 
accuracy of the numerical method. If the numerical solution converges, the relation between 
the different numerical solutions and the real one will be: 

F real = F h + 0(h n ), (38) 

where n is the order of convergence and F h is the numerical solution with step size h. 
For halved step sizes the error is reduced correspondingly: F rea i = F h l 2 + 0((|) n ) and 
Freai = F h / 4 + 0((f ) n ). By defining (as in [18]): cl = F h - F h ' 2 and c2 = F h l 2 - F h /\ we 
find the convergence rate: 

n = log 2 (^j ■ (39) 

Figure 4a depicts n. The convergence rate is approximately linear n ~ 1 or higher for 
almost the whole domain. However, n diverges at some points. This arises from crossing 
of s = const, lines as can be shown in figure 4b. The upper panel of this figure shows a 
one dimensional projection of n along constant u-ray: u = 0.4, and the lower panel shows 
the corresponding field function F = s(u, v) along the same constant u ray for different grid 
densities. 
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(a) The convergence rate - n. (b) One dimentional convergenge along a constant u-ray, 

u = 0.4. The convergence rate - n (upper panel) and the filed 
function - s for different grid densities (lower panel) 

FIG. 4: Convergence for a slightly subcritical run with a = 10~ 4 . The convergence rate - n - in panel (a) 
is derived from the field function s(u, v), with 2 16 , 2 17 and 2 18 grid points in u and v directions (blue, red 
and green lines correspondingly in s plot on panel (b)). In panel (b) the divergence of n around v ~ 0.46 
and v ~ 0.63 is caused by s lines crossing, as shown in the zoom window. 

IV. RESULTS 

The first feature of the collapsing field is the formation, or not, of a black hole. To 
examine this we plot a diagram of v vs r, the area coordinate, for different values of u (See 
figure 5). Each line represents an outgoing null ray of a constant u, namely a trajectory of 
a photon emitted from the origin at u. Rays with u small enough show a flat-like spacetime 
behavior, for which v ~ u + 2r. For small values of u, at early times, the rays don't encounter 
a strong gravitational field and they escape to infinity. At later times the gravitational field 
becomes stronger and the outgoing rays are bend more and more before they eventually 
manage to escape. Once a black hole forms these rays are trapped. For subcritical initial 
configuration the field disperses and all outgoing null rays reach infinity (see figure 5a). For 
a supercritical initial configuration a horizon appears when an outgoing null ray doesn't 
escape and doesn't reach future null infinity but rather remains in the same radius r for all 
values of v. Later rays that emerge from the origin collapse back to the origin (see figure 
5b). 
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(a) Subcritical collapse, p < p*. (b) Supercritical collapse, p > p*. 

FIG. 5: Outgoing null rays: v vs r. Both plots have the same parameters except for the initial field 
amplitude: a = 10~ 5 , v rnax — 0.768, the grid density is 2 15 points in each direction. 

While, as expected the Gauss-Bonnet terms don't change the overall classical GR behav- 
ior, the critical behavior is lost. For classical GR in D = 5 dimensions the black hole radius 
scales as R oc (p — p*) 1 . Figure 6 depicts the black hole radius as a function of \p — p*\ for 
supercritical evolutions with a = 10~ 4 . Note that we use here the scaling relation of the 
Schwarzschild radius instead of the mass (see section II F). The dependence of the black hole 
radius on the initial amplitude is monotonic, i.e. for larger values of p the black hole radius 
is larger. However in the EGB gravity we don't observe the classical power-law relation. 
The existence of the scaling relation is related to the self similar properties of the classical 
solution [2, 19]. Since self similarity is not preserved in EGB gravity we expect that the 
scaling relation will also be violated. 

In classical GR the solution is discretely self similar (DSS) just below the black hole 
threshold. The field reaches the origin, oscillates and then disperses or collapses, depending 
on whether it is subcritical or supercritical. These oscillations don't depend on initial con- 
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FIG. 6: The black hole's Schwarzschild radius - r s vs. the difference between the amplitude and the 
critical one - \p — for supercritical cases with a = 10~ 4 . In the classical GR solution this relation is a 
power law. Here, in EGB gravity, a power law is not observed. 



ditions and they decay in a DSS pattern (see figure 2a). The addition of the Gauss-Bonnet 
terms destroys this behavior. With these terms the field still oscillates, but the self similarity 
disappears (See figure 8). The oscillations grows at first as if the additional terms in the 
Lagrangian amplify the field, prevent it from decaying and keep the oscillations alive for a 
longer time. Eventually the oscillations decay and the field disperses. 

In the classical GR solution, we observe more and more self similar transients when 
approaching the critical amplitude from both sides. In EGB solution, while approaching the 
critical amplitude more and more oscillations are also observed. However, these oscillations 
are not self similar and their scale does not decrease. Figure 8 depicts the contours of the 
scalar field s for a set of subcritical solutions (panels (a-d)), with a = 10 -5 , with growing 
amplitudes approaching the critical one, and a set of supercritical solutions (panels (e-h)) 
with amplitudes decreasing towards the critical one. An increasing number of oscillations is 
observed as p approaches p*. The inserts in the supercritical solution, figures 8 (e-h), depict 
the v — r diagram, demonstrating black holes formation. A cutoff in field diagram is a sign 
for a singularity (see the scheme in figure 1). 
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FIG. 7: The typical wavelength of scalar field oscillations in EGB gravity near the black hole threshold 
vs. the coupling constant a. Red line indicates a linear fit with a slope: m — 0.54 ± 0.05. Namely, the 
wavelength is proportional to \fa as expected from the dimensional analysis. 



Obviously the self similarity disappear not only in the field functions, but in all metric 
functions and their derivatives. Particularly interesting is the behavior of the Ricci scalar 
(equation (31)) as seen in Figure 9 panels (a3) and (b3), and Figure 10 panels (c3) and (d3). 
This figures show the Ricci scalar at the origin, Rir = 0), as a function of u for different 
values of a. Purple color indicates negative values and red indicates positive values of the 
Ricci scalar. Figure 9(a3) shows the classical solution, i.e. a = 0. In classical GR the 
condition a < (equation 32) is always satisfied, thus R(r = 0) is always negative and it 
never changes sign. On the other hand for EGB gravity this condition is inevitably violated 
for amplitudes close enough to the critical one (see figures 9(b3), 10(c3) and 10(d3)). As p 
approaches p*, the metric function a tends to zero and the derivative of the field z grows. 
At some point the condition (32) is violated and the Ricci scalar changes sign, indicating 
a change in the local geometry. This demonstrates that near the black hole threshold for 
a > the local geometry is different from the classical GR geometry. 

Figures 9 and 10 compare slightly subcritical solutions with different values of the coupling 
constant, a = 0, 10~ 6 , 10~ 5 , 10" 4 in each column (columns a,b,c and d). The first, upper row 
of each figure (panels al, bl, cl and dl) shows contours of the field function s. The second 
and the third rows (panels a2-3, b2-3, c2-3 and d2-3) show the field function s and the Ricci 
scalar on the origin (r = 0) vs. u. The 4 th row (panels b4, c4 and d4) shows a contour 
plot of \aR\. It indicates the regions where the additional curvature terms are significant, 

18 




FIG. 8: Contours of the scalar field function s in EGB gravity with a coupling constant a = 10~ 5 , grid 
density 2 15 , v max = 0.768, v c = 0.22. The initial field amplitude increases or decreases towards the critical 

amplitude: p — > p*. Plots (a-d) - subcritical collapse: p < p*. Plots (e-h) - supercritical collapse p > p*. 
The inserts demonstrate the formation of black holes as seen in the v — r diagram. The initial values of the 
field amplitude are: p = 0.1654(a), 0.16557(b), 0.1657(c), 0.165723(d), 0.1676(e), 0.166(f), 0.1658(g), 

0.16573(h). 

i.e. regions where |<xR| > 1. Naturally with larger values of a this region grows and \aR\ 
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reaches larger values. 

The plots of the field function - s - on the origin (b2, c2 and d2) nicely show a "beat "-like 
pattern, increasing and then decreasing, in the field pulsation's strength. However, for small 
values of a, (see figure 9(b2)), the field behavior on the origin resembles, at least initially, 
a self similar behavior in regular GR. This could be explained by the low and insignificant 
values of the higher order terms in these regions (see figure 9(b4)). At the same time a 
comparison of the field s and the Ricci scalar at the origin for different values of a (panels 
(a2-3), (b2-3), (c2-3) and (d2-3)) reveals that the new "beat" form of the field oscillations 
appears at the same retarded time u at which Ricci curvature changes sign and becomes 
positive. 

The scalar field oscillations (figures 8, 9 and 10) show a typical wavelength which depends 
on the value of a (see figure 7). All the lengths are measured in the simulation length units 
[u\. As expected from a dimensional analysis the typical wavelength of the oscillations is 
proportional to y/a. 

V. SUMMARY 

We have developed a numerical scheme for simulating the dynamical collapse of a spheri- 
cally symmetric massless scalar field in EGB gravity. This model for gravity includes higher 
(quadratic) order curvature corrections to the Hilbert action. These corrections induce 
changes in Einstein equations, which govern the evolution of the system when the curvature 
is large. 

We find that the addition of higher order curvature correction destroys the classical Chop- 
tuik phenomenon. The introduction of the dimensional coupling constant a, which has a 
units of length 2 , destroys the scale invariance of the system. As a consequence the self 
similar behavior, which is an integral part of the critical phenomena in regular GR, disap- 
pears. Instead the solution shows a different pattern of pulsations with a typical wavelength, 
which is proportional to y^, as expected from a dimensional analysis. The changes in the 
oscillations pattern are accompanied by changes in the sign of the Ricci scalar at the origin, 
indicating a change in the local geometry of the space-time. 

We thank Shahar Hod and Nathalie Deruelle for many helpful discussions, Stanley Deser 
and Hideko Maeda for useful remarks and Yonatan Oren and Evgeny Sorkin for assistance 
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with the numerical calculations. 




FIG. 9: A slightly subcritical collapse with different values of the coupling constant a. Panels (al-3) 
correspond to a classical GR solution, a = 0. A self similar behavior can be observed. Panels (bl-4) 
correspond to EGB gravity with a = 1CP 6 . The first row, panels (al) and (bl), presents contour plots of 
the field function s. The second row, panels (a2) and (b2), shows the field function s at the origin (r = 0) 
vs. u. The third row, panels (a3) and (b3), shows the Ricci scalar on the origin vs. u in a logarithmic 
scale. A red color indicates positive values and a purple color negative values of R. Panel (b4), displays 
contour plot of |<xR|, showing the regions where the higher order terms are significant, i.e. \aR\ > 1. 
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